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ABSTRACT 

While supermassive black holes (SMBHs) play an important role in galaxy and cluster evolution, 
at present they can only be included in large-scale cosmological simulation via subgrid techniques. 
However, these subgrid models have not been studied in a systematic fashion. Using a newly developed 
fast, parallel spherical overdensity halo finder built into the simulation code FLASH, we perform a 
suite of dark matter-only cosmological simulations to study the effects of subgrid model choice on 
relations between SMBH mass and dark matter halo mass and velocity dispersion. We examine three 
aspects of SMBH subgrid models: the choice of initial black hole seed mass, the test for merging two 
black holes, and the frequency of applying the subgrid model. We also examine the role that merging 
can play in determining the relations, ignoring the complicating effects of SMBH-driven accretion and 
feedback. We find that the choice of subgrid model can dramatically affect the black hole merger rate, 
the cosmic SMBH mass density, and the low-redshift relations to halo properties. We also find that 
it is possible to reproduce observations of the low-redshift relations without accretion and feedback, 
depending on the choice of subgrid model. 

Subject headings: black hole physics, cosmology:theory, dark matter, galaxies: evolution, large-scale 
structure of universe, methods: numerical 



1. INTRODUCTION 

Supermassive black holes (SMBHs) play a number of 
roles in the evolution and dynamics of galaxies and clus- 
ters of galaxies. These black holes, with masses of at 
least a million solar masses, and their a ssociated accre- 
tion disks drive quasars at high re dshift (j Fan 2006 ?) , reg - 
ulate star formation in galaxies ([Hopkins et al.l 12010) . 
inject therma l energy into t he intracluster medium via 
powerful jets (|Gu et al.ll2009f ) . and may even play a large 
role in establishi ng kiloparsec-scale micr ogauss magnetic 
fields in clusters ([Carilli fc Tavlorll2003) . 

Of particular interest arc the correlations discov- 
ered between SMBH mass and other observed quan- 
tities of galaxies. The first discovered relation- 
ship , between SMBH mass and bulge stellar luminos- 
ity ([Magorrian et al.l [19981. was intriguing but suffered 
from large scatter. Subsequent searches found a tight 
correlation between SMBH mass and bulge veloc i ty dis - 
persion a (|Tremaine et al.l l2002t iGultekin et al.l 120091) . 
although comparable improvements have also been mad e 
in the correlations to stellar luminosity ([Graham! I2007j ). 
Most recently, black hole mass has been linked to dark 
matter halo mass M to t, both indire ctly by measur - 
ing galactic circu lar velocities (e.g., ([Ferraresel [2002); 
(jBaes et al.ll2003D ) and b y direct estimates of halo mass 
via gravitational lensing ([Bandara et al.l 12009). The lat- 
ter measurements agree well with theoretical predictions 
of the relationshi p between ha lo virial mass and galactic 
circular velocity (jCroton 20091) . 

These relationships imply a correlation between the 



growth of SMBHs and that of their host galaxies. Since 
structures form in the universe via h ierarchical cluster- 
ing of smaller objects (|Baughi feOOSf ) , black holes car- 
ried along with their ho sts should tend to merge as 
well ([Hopkins et al.l 120051 ). This may provide a simple 
and direct scaling betw een halo and black ho le mass, es- 
pecially at low redshift ([Volonteri et al.ll2~003l ) . However, 
feedback processes driven by accretion disk systems may 
also contribute to the observed correlations: greater ac- 
cretion can lead to larger feedback events, thereby reduc- 
ing the accretion rate and coupling the black hole ma ss 
to the surrounding system ([Cattaneo fc Tevssierll2007t ). 

Attempts to explain the observed relations via cos- 
mological simulations require large volumes (to gather 
enough objects) and high resolution (to capture galaxy- 
sized structures). However, even with current computing 
resources, large- volume simulations are unable to cap- 
ture all the intricate physical processes that dominate 
in the formation, merging, and feedback of black holes. 
Hence, these processes must be added in post-pro cessing 
as semi-analytic models (e.g., iMicic et al.ll2007f) or in- 
clude d in-situ as subgrid models (e.g.. iBooth fc Schavd 
120091 ; iSiiacki et all 120071 ; iDi Matteo et al.ll2008l ). " 

Many processes have been proposed to explain the for- 
mation of seed SMBHs in the early u niverse, includ- 
ing remnants of Popula tion III stars ([Madau fc Reesl 
120011; iWise fc Abel 120051). direct collapse of gas in cen- 
tral bulges ([Koush iappas et al.l l2004t [Bcgclma n et al 



2006), and merging of smaller black holes ([Islam et a! 
I2004D . Merging black holes are difficult systems to 
study, si nce they interact with their gaseous 
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ronment (Mayer et alj l2007f) . emit gravitational rad: 
ation (Scsan a et a 1.1 I2004D . and can suffer kicks due 
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to merging (|Baker et all |2008). Additionally, the 
self-regulating feedback processes emerging from ac- 
cretion onto black holes and the subsequent for- 
matio n of jets and bubbles a re not fully under- 
stood (jVernaleo fc R eynolds 200(3) , es pecially when con- 
sider ing the effects of m agnetic fields ([ Ruszkowsk i et al.1 
120071) and turbulence (jBruggen fe Scannapiecol 12009). 
Consequently, subgrid models must make many simpli- 
fying assumptions in treating these processes. 

The variety of plausible scenarios for forming and 
merging black holes and applying feedback processes al- 
lows modelers great latitude in developing and adjusting 
models to fit observations. Universally, all aspects of the 
formation and evolution of SMBHs are combined in the 
same simulation. However, we believe that subgrid mod- 
els of the initial seeding and merging of SMBHs (which 
are linked in subgrid models to the properties of the sur- 
rounding dark matter) should be separated from models 
of accretion and feedback (which depend on the local gas 
physics). This way, we can better understand the role 
that merging alone plays in developing the Mbh — cr and 
Mbh — A/tot relations and the effects of changing subgrid 
models on those same relations. 

Thus, in this paper we examine dark matter-only sim- 
ulations of the growth of structure in a cosmological vol- 
ume including subgrid models to track the formation and 
merging of SMBHs. By comparing the results of several 
plausible scenarios for models against observed relations, 
we will determine how much of those relations is due 
to low-redshift and large-scale evolution of SMBHs, and 
how these models may affect the final outcomes, inde- 
pendent of any gas accretion or feedback. 

To follow our dark matter halos, we have developed 
a new, parallel, fast halo finder built directly into the 
simulation code FLASH v2.5. FLASH is an adaptive- 
mesh re finement (AMR) cod e for astrophysics and cos- 
mology (|Frvxell et al.l I2000D . FLASH solves the N- 
body potential problem with a particle-mesh multigrid 
fast Fourier transform method (Ricker 2008). It uses 
smoothed cloud-in-cell mapping (Ricker et al. in prepa- 
ration 2010) for interpolating betwe en the mesh and par- 
ticles (Hockncy & Eastwood 119881 ) and a second-order 
leapfrog integration scheme for variable time step parti- 
cle advancement. 

In the following section we discuss the precision and 
valid mass ranges for our new halo finder. In Section [3] 
we outline the numerical aspects of our approach and 
the black hole formation and merging subgrid models 
employed. Finally in Section f?] we compare our results 
to observations of the Mbh— a and Mbh.— -Wtot relations to 
test the validity of the models. Additionally, we provide 
a discussion and analysis of the performance and parallel 
scalability of our halo finder in the Appendix. 

2. THE HALO FINDING METHOD 

We base the halo finder in FLASH on a spherical- 
overdensity (SO) technique. Throughout, we will identify 
our new halo finder by "pSO" , for parallel spherical over- 
density. In this approach, halos are defined by spherical 
regions within which the mean density is greater than 
some defined threshold. We begin by mapping parti- 
cles onto the simulation mesh with smoothed cloud-in- 
cell mapping and identifying peaks by finding zones with 
densities greater then all surrounding zones and greater 



Table 1 

Parameters controlling the pSO halo finder in FLASH. 



Parameter 


Description 


Value 


Ap ca k 


Overdensity for identifying a halo center 


200p crit 


^search 


Overdensity for defining a halo 


200p crit 




Criterion for completing a radius search 


0.2Ax 


Airmail 


Criterion for aborting a radius search 


0.5Az 




Minimum resolvable halo radius 


l.OAx 


Note. - 


Ax is the uniform grid resolution. 




than A pea 


kPcrit- Here and throughout, p cr j t refers to the 



comoving critical density of the universe, 
3H 2 

Pcrit = [%+^(l + ^ 3 ] , (1) 

where subscripts of here and throughout refer to 
present-day values. The zone midpoints serve as centers 
of potential halos. Using a binary search, we compare 
the average density within the current search radius to 
Asearchpcritj selecting new search radii appropriately. Ini- 
tially, search radii are doubled until the enclosed density 
is below the threshold. Only then does the binary pro- 
cedure begin. When two successive search radii differ by 
no more than a chosen amount, defined by the parame- 
ter Ai? st0 p, we stop the search. If, during the search, a 
radius falls below a cutoff value, Ai? sma n, we the abort 
the search and disregard the halo. Finally, we remove 
any finished halos that have radii smaller than Ai? m j n 
from the catalog. We also remove satellite halos whose 
centers fall within the radii of larger companions. This 
leads to a more consist ent mass function and conserves 
halo mass (White 2002) . We may also optionally remove 
satellite halos which intersect larger neighbors. While 
this must be done to strictly conserve halo mass, many 
authors include these satellite halos, since it is sometimes 
us eful to identify satellite structures (see the comparison 
in lEvrard et all (|2008f ) for a discussion of this decision). 
Table [1] lists the parameters controlling our halo finder 
and our chosen values. In the table and throughout, 
Ax is the grid resolution. Note that for simulations in 
which the halos may span an adaptively refined region, 
Ax will refer to the highest-resolution uniformly refined 
mesh level. 

For this work, we chose A pca k = A scarc h = 200. How- 
ever, other values may be chosen for other uses of the 
halo finder. For example, A soarc h = 500 would be appro- 
priate for generating mock o bservations of X-ray cores for 
comparison to observations (jEvrard et al.|[l996D , while a 
smaller value might be useful for triggering refinement. 

To evaluate our halo finder, we performed dark matter- 
only simulations in a cubic 128 /i -1 Mpc box with 256 3 , 
512 3 , and 1024 3 zones, giving resolutions of Ax = 500, 
250, and 125 ft,~ 1 kpc respectively. All runs used 256 3 
particles, giving a mass resolution of 1.3 x 10 10 M Q . 
For these tests, we chose cosmological parameter val- 
ues of Q M ,o = 0.26, fi A ,o = 0.74, and H = lOO/i = 
71 km s~ 1 Mpc~ 1 for comparison to runs used in previous 
works. We compared our pSO halos against spherically 
overdense regions drawn from a friends-of-fr iends (FOF) 
halo finder with a linking length of 0.2 (see iLukic et al.l 
120091 for an analysis of SO halos drawn from an FOF 
catalog). The halos drawn from the FOF catalog used 
the most-linked particle as the halo center. Halo radii 
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Figure 2. Error in matched halos before and after recentering 
FOF halos on pSO centers, as a function of halo radius. The error 
is defined in the text. 



Figure 1. Demonstration of different choices of halo center due to 
recent merging. The thin lines indicate projected density contours 
of the FOF halo. The interior contour lines define thresholds of 
50 and 200 particles per zone, while the outermost contour shows 
a threshold of 1 particle per zone (i.e., the boundary of the FOF 
halo). The two thick circles show the boundaries of the spheri- 
cally overdense regions identified by pSO. The separation of the 
larger pSO halo center (circle) and the FOF halo center (square) 
is roughly one half a zone. 

in the FOF catalog were determined by starting with a 
large radius and moving inwards at very small increments 
(much smaller than the Ai? stop used in pSO) until the 
interior density exceeded the same threshold as used by 
the pSO halo finder. This approach is more precise, but 
much slower, than the binary search technique used in 
pSO. We use a linear search here since with such high 
precision a binary technique may not complete. 

We thus have several possible sources of differences be- 
tween the SO and FOF halo catalogs: (1) pSO chooses 
the zone midpoints as the halo centers, while FOF halos 
use the most-linked particles, (2) the incremental search 
radii when finding spherical overdense regions in the FOF 
catalog are much smaller than those used in pSO, (3) 
FOF tends to find more small halos than the grid-based 
peak finding in pSO, since smaller halos may be irregu- 
larly shaped, and (4) FOF will tend to bridge two nearby 
halos, even if they have distinct spherically overdense re- 
gions. 

Figure [T] demonstrates both the problem of selecting 
a halo center and that of pSO counting more satellite 
halos than FOF. Shown is a ~ 10 14 /i _1 M Q halo drawn 
from the run with 250 /i -1 kpc resolution. This halo is 
undergoing a merger with a smaller satellite, and thus 
FOF is bridging two distinct spherical regi ons into a sin- 
gle ha lo. In agreement with the findings of lEvrard et al.l 
( 2008), we find roughly ~ 15% of our halos are satel- 
lites. While the most-linked particle approximates the 
potential minimum well, the pSO halo center, which is 
simply the zone midpoint, is within a zone radius of this 
point. The bridging effect leaves the smaller pSO halo 
unmatched. 

We match halos in the pSO and FOF catalogs by find- 
ing intersections. Any halos from the two catalogs that 
overlap such that they contain each other's centers are 



considered a match. Due to the ambiguity of halos with 
very few particles, we only compare halos with more than 
100 particles. With this procedure alone, about 15% of 
the pSO halos remained unmatched. To determine if this 
was due to the FOF bridging effect, we removed from 
consideration any pSO halo that was within a linking 
length of an already matched FOF halo. This accounted 
for all the unmatched pSO halos. 

To account for the effects of choice of halo center, we 
re-computed overdensities in the FOF catalog using the 
potential minimum halo centers found in the matching 
pSO catalog. Note that we found that the differences be- 
tween halo centers were always smaller than a zone spac- 
ing. There may be some cases, however, where a large 
asymmetry will separate the most-linked particle and the 
maximum density peak further than a zone spacing, es- 
pecially at very high resolutions. Here, asymmetries will 
manifest themselves in a higher relative error between 
matched halos. Figure [2] shows errors in the matched 
halos in the 125 h kpc catalog at z — before and af- 
ter recentering. We define the errors as the difference in 
mass divided by the sum of the la uncertainties in the 
halo mass: 



We estimate the halo mass uncertainty by assuming an 
uncertainty in the halo radius of 0.5Ax. Thus, assuming 
constant halo density, the mass uncertainty is 

<tm = -— AxcxM 2/ 3 , (3) 
Z K 

where M and R refer to the halo mass and radius, 
respectively. For halos above 10 14 M Q in mass, this 
le ads to oyilM k, 0.1, which agrees with the estimates 
of Bhattachar va et al.l (|20100 . While most matched ha- 
los have small error, a few differ by as much as 0.7, es- 
pecially at lower masses. However, most of this is due to 
the bridging effect's having moved the most-linked parti- 
cle away from the potential minimum. After recentering, 
errors for all halos larger than a zone radius drop to be- 
low 0.1. The small gap in halo sizes near i?20o/-Rzonc = 
is due to our choice of Ai? st0 p and the binary search pro- 
cedure. 
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Figure 3. Error in matched halos at different simulation resolu- 
tions. Bold lines denote average error, while thin lines describe 
maximum error for each mass bin. The error is defined in the text, 
and is computed after recentering. 

Also from Figure [2] we find that even after recentering, 
halos smaller than a zone do not match well to FOF 
halos. Thus we set the parameter Ai? m i n to l.QAx and 
reject any halo smaller than this. Given a resolution, 
this sets our minimum resolvable halo mass. Note that 
this criterion is consistent with the estimate f or minimum 
resolv able mass for FOF halos as defined bv lLukic et al.l 
(pOOl : 

n min = A ^ m7 ! P ! n i 3 PW1 + z f + fiA,o] » (4) 

where n g and n p are the number of zones and particles 
per side, respectively, and A = 200 is our chosen over- 
density. 

Figure [3] shows maximum and average errors after re- 
centering for each of 10 mass bins at z — for each 
simulation resolution. With our approach, we are able 
to maintain average errors of less than 0.04 for all resolv- 
able halos, while a small number (< 5%) of halos have 
maximum errors of up to 0.12. Although the average er- 
ror stays consistent across all masses, the maximum error 
varies as much as 0.04 for adjacent bins, and it rises with 
decreasing mass. This behavior is due to a small num- 
ber of irregularly-shaped halos. Also, since smaller ha- 
los have smaller uncertainties, the errors tend to become 
larger with smaller halo mass. 

Figure |4] shows errors after recentering at various red- 
shifts in the Ax = 125 /i _1 kpc run. Across all resolvable 
masses and redshifts from z — to 2, we are able to 
maintain average errors less than 0.05. However, max- 
imum errors in the smallest mass bins reach as high as 
0.18. Again, the variability of the maximum errors is due 
to a small number of halos. 

By adjusting Ai? stop , we are able to reduce both the 
average and maximum errors. However, at smaller values 
we found that the binary search procedure had difficulty 
converging on a value for some halos, and the halo finder 
ran for an unacceptable amount of time. Larger values 
produced unacceptably high errors. For the value we 
chose, after accounting for different choices of halo cen- 
ter, the halos produced by our new pSO halo finder are 
statistically indistinguishable from matched halos drawn 
from a traditional FOF halo finder. 



Figure 4. Error in matched halos at different redshifts for run 
with 125 h _1 kpc resolution. Bold lines show average error for each 
bin, while thin lines denote maximum error. The error is defined 
in the text, and is computed after recentering. 
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Figure 5. Mass functions at different resolutions compared 
against halos drawn from the 125 h ~ kpc resolution FO F catalog 
and against the mass function from Warren et al. ( 2 009 ). 

Figure [5] shows the mass function for the pSO halo 
catalogs at the three resolutions compared to SO ha- 
los drawn from an FOF catalog with a resolution of 
125 /t -1 kpc and c o mpar ed to the mass function obtained 
by iWarren et al.l (|2006t ). As expected, pSO captures 
more halos as the resolution increases. However, pSO 
captures fewer halos near the resolvability limit than 
FOF does. The procedure for mapping and smoothing 
tends to lower the central density, especially with halos 
near the resolvability limit. Thus our criterion for select- 
ing a peak, Ai? pca k = 200, may be too stringent. We 
found that lowering this value captures more halos, but 
at the expense of identifying too many potential halos 
that end up below the resolvability limit. Even at values 
as low as 100, we still could not replicate the FOF mass 
function at these masses, while adding roughly four times 
as many candidate halos that were ultimately rejected. 

Finally, to speed up processing we may increase 
Ai? sma n. However, the binary search may briefly fall 
below this threshold before settling on a larger, and cor- 
rect, radius. Thus we want to pick the largest possible 
value that does not cause us to remove resolved halos. 
We found that at any values above 0.5Aa; we began to 
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reject valid halos before their search had completed. 



3. SUBGRID MODELS 

We identified three areas in which subgrid models may 
differ: the initial mass of seed black holes, the prescrip- 
tion for determining if two black holes merge, and the fre- 
quency of finding halos, creating black holes, and check- 
ing for mergers. For each of these aspects we study two 
possibilities: a "pessimistic" and an "optimistic" sce- 
nario. This leads to a total of eight combinations of 
models. 

For the initial mass of seed black holes, au- 
thors who have performed cal cu lations similar to 
ours (e.g.. iBooth fc Schavd 120091 ; iSiiacki et all 120071 
; iDi Matteo et alll2008l ) typically choose a constant seed 
mass of 10 5 M Q , and this is the value we will choose for 
this model. The value chosen is typically motivated by 
a desire not to vio late constraints on th e observed black 
hole mass density (jShankar et all 120041) , especially when 
black holes are allowed to accrete a significant portion of 
their final mass. 

However, this approach may have several weaknesses. 
First, black holes with masses as high as 10 9 M^ are in - 
ferred to exist from quasar activity at z ~ 6 (|Fanl f2006). 
Second, depending on the frequency of halo finding, when 
halos are initially detected they may not have identical 
masses, and if whatever physical processes which cre- 
ate the observed relations are already present, then the 
seed black holes should scale with halo mass. Finally, at 
lower redshifts the halo finding algorithm may spuriously 
merge two halos. If we instantly merge black holes and 
in the following step the halo finder identifies two sepa- 
rate halos, one will be without its black hole. It may be 
inappropriate to re-seed the halo with a low-mass black 
hole. 

Alternatively, we may derive the seed black hole mass 
from a form of the observed Mbh — o relation. This ap- 
proach may seem tautological; however, accretion and 
feedback processes may dominate at hi gh redshift and 
at scales below wh ich we can resolve (Shankar 2009; 
iMerloni et alj |2010). The observed Mbh — v relation at 
z = may then simply be a consequence of pure merg- 
ing at late times and large scales. Also, this approach 
may be useful in testing other SMBH-related processes 
and relations in groups and clusters. However, this ap- 
proach is complicated by the fact that calculations of a 
directly in simulations suffer from high scatter, especially 
at high redshift when halos first form and at low redshift 
when satellite halos may be contaminated by hot parti- 
cles belonging to a larger neighbor. Fortunately, we may 
seed halos by relating the black hole mass directly to the 
halo mass through the observed relation (Bandar a et all 
2009): 



"*(^) =(8.^°.ll) + (1.55 ± 0.31)l„ s ( I5 ^ 

This relation comes from observations at z = 0.2 of bulge 
velocity dispersions and halo masses using gravitational 
lensing. These measurements are linke d to black hole 
mass by assuming the Mbh — f relation of lGultekin et all 
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Figure 6. ctdm as a fun ction of halo mass a t redshift z = 0.2, 
compared to observations jBandara ct al. 2009) and best fits from 
simulations (Evrard ct al. 2008). The latter is an extrapolation 
below ~ 10 14 M G 
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where ctq = 200 km s . We will assume the relation in 
Eq.© holds to high redshift. 

Figure |B] shows velocity dispersion as measured in one 
of our runs (detailed below) as a function of mass at 
z = 0.2. While we agree with observations and the best 
fit fro m a suite of simulations described by lEvrard et all 
( 2008), there is significant scatter, especially at low mass. 
Note that since we are performing dark matter-only sim- 
ulations, we do not directly measure the velocity disper- 
sion of the galactic central bulge; instead we measure 
that of the dark matter in the entire halo, defined by: 
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where a(t) is the scale factor, N p is the number of parti- 
cles in the halo, Vij is the jth velocity component of the 
zth particle, and Vj is the jth component of the center-of- 
mass velocity. Note that velocities in the simulation are 
comoving q uantities. Thi s plot and other results (as sum- 
marized by ICrotonl 120091 ) indicate that cdm and ab u i ge 
agree to within a factor of order unity, so we may sub- 
stitute one for the other. 

Next, we may change the prescription for testing black 
hole mergers. Most approaches to studying the growth 
of SMBHs r ely on halo merg er histories derived from 
sim ulations (iMicic et al.1 12007|) or Press-Sch echter mod- 
els (|Menou et al.ll2001t iVolonteri et alll2008l ). These ap- 
proaches assume that each halo contains a single SMBH, 
and whenever two halos merge their respective black 
holes instantly merge as well. While this approach is 
highly optimistic, it does provide an upper limit to pre- 
dicted merger rates, an important constraint on upcom- 
ing g ravitational wave experiments (|Amaro-Seoane et all 
2007). It is also easy to implement and provides a simple 
way to compare results from simulations to approaches 
based on merger trees. However, this approach can lead 
to spurious black hole mergers, since halos may tern- 
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porarily intersect without merging. We will include this 
scenario in our study. 

Subgrid models us ed directly in simulations (e.g., 
iBooth fe Schavd 120091 ) have used a more sophisticated 
approach in merging black holes. In these simulations, 
SMBHs only merge when they are within some defined 
distance (for Lagrangian codes, this is typically the soft- 
ening length) and have relative velocities below some 
threshold. This is to avoid merging black holes that are 
only passing by each other and are not part of a true co- 
alescing system. The velocity threshold varies by author, 
but it is usually taken to be the local gas sound speed 
or the circular velocity near the larger black hole of the 
merging pair. However, we found this test to be overly 
restrictive; without gas to slow down black holes in clus- 
ter cores, merging times inferred from these criteria are 
larger than the Hubble time and hence almost no black 
holes merge using this model. Since observations indi- 
cate that black holes do merge (Mcrritt & Milosavlievic 
2005j), we must apply a less stringent test. 

Simulations of merging galaxies suggest that the time 
for black holes to mov e from kiloparsec to parsec scales 
is typically ~ 10 Myr (jDotti et al.ll2007t ). T his assump- 
tion w as used in the merger-tree analysis of iMicic et al.l 
(2008). When applying thi s test we will ignore th e ef- 
fects of gravitational recoil (Bogdan ovic et aT1l2007f) and 
the "final parsec problem" ([Berczik et al.ll2006f) . Thus, 
in this merging test we will only merge black holes if 
they are both within the same halo and are within two 
grid zones of each other, and their relative velocity times 
10 Myr is less than that same distance. This approach 
allows black holes to pass near each other without merg- 
ing if they are not part of a truly merging system and 
accounts for the time needed for the black holes to merge 
below our resolvable scales. 

Finally, we may alter the frequency of performing our 
subgrid analysis. One common approach is to check for 
new black holes and allow mergers such that the inter- 
val between successive c hecks is evenly spaced in the log 
of the expansion factor. IBooth fc Schavd ()2009f ) employ 
the shortest such interval, such that a nex t = l-02a current . 
This approach requires roughly 1 /6 the number of halo 
searches compared to searching every time step. Since 
our halo finder is designed to be inexpensive, we may 
perform checks at every time step. This may cause spu- 
rious formation and merging of black holes, especially 
with small halos near the resolvability limit. However, 
with infrequent checks we may miss the formation and 
merging of new small halos, underestimating both the 
amount of black hole mass and the merger rate in the 
simulation. We will study checking both every time step 
and at an interval of A log a = log 1.02. 

Table [2] summarizes the aspects of the subgrid model 
we are studying, our choices for modifying each aspect, 
and a shortened name that we will use to identify the 
models in plots and tables. For example, a model that 
uses a constant initial seed mass, merges black holes in- 
stantly, and performs a check at every time step would 
be identified by "con, halo, dt" . We note that the most 
commonly used model in the literature uses a combina- 
tion of uniform initial mass, velocity tests at every time 
step for merging, and seeding new black holes evenly in 
log expansion factor. For our study, we have combined 
merging tests and seeding in the same step. This is closer 



Table 2 

Aspects of Subgrid Models Used in the Creation and Merging of 
SMBHs. 



Aspect 


Model 


Short Name 


Seed mass 


Constant 


con 




A^bh — c relation 


m-s 


Merging strategy 


Instantly on halo merger 


halo 




Distance and velocity test 


prox 


Frequency 


Every time step 


dt 




Evenly spaced in log(a) 


log 



to the approach used in merger tree analysis, where anal- 
ysis can only take place on the available halo catalogs. 



4. COMPARISON OF MODELS 

For all calculations, we used concordance parameter 
values of M ,o = 0.238, Q A)0 = 0.762, and H = lOO/i = 
73.0 km s" 1 Mpc -1 . All runs took place in a three- 
dimensional box measuring 50 ft, Mpc per side with 
512 3 particles and 1024 3 zones per side, giving a mass 
resolution of 6.15 x 10 7 h Mq and spatial resolution 
of 48.8 ft, _1 kpc. There was no refinement of grid spac- 
ing. All simulations used the same initial conditions: 
unperturbed particle positions were situated on a grid, 
and the initial velocities and positions were perturbed us- 
ing Gaussian fluctuations normalized to cr§ = 0.74. We 
assumed P(k) fr om a ACDM cosmo logy. We used the 
GRAFIC2 code (jBertschingeri I2001D to generate these 
initial conditions. All computations started at a redshift 
of z = 56.8. 

Using our halo finder, our minimum resolvable halo 
mass is ~ 10 10 M©. To help avoid spurious mergers, 
we did not include satellite halos (halos that intersect a 
larger neighbor) in the halo catalog used by the seed- 
ing and merging models. We seed black holes in any 
resolvable halo by creating a black hole particle in the 
simulation with a dynamical mass equal to the mass of 
the black hole and velocity equal to the center-of-mass 
velocity of the parent halo. When we merge black holes, 
we remove the smaller of the pair and add its mass to 
the larger of the pair. The larger SMBH maintains its 
position and velocity. We do not re-center black holes on 
halo potential minima. 

While re-positioning black holes at potential minima 
reduces ejections due to scattering, this process may ar- 
tificially promote merger rates. With our SMBH setup, 
we do indeed see a large ejection rate. Fortunately this 
is not an issue for our analysis: the vast majority of ejec- 
tions are of un-merged black holes from low-mass halos 
at high redshift. Consequently, small black holes are sim- 
ply replaced by another small black hole, and the merger 
rates and relations remain unaffected. Obviously, this 
would pose a problem for models that include accretion. 
The ideal solution in this case would be to initialize on 
halos well above the minimum resolvable halo mass and 
to adequately smooth the gravitational potential near 
the SMBH so that it does not experience significant two- 
body effects. Also, by using the black hole mass as the 
dynamical mass, we may underestimate the dynamical 
friction. However, we have found that almost all black 
holes lie within a zone of their host's potential minimum, 
so this does not affect the merger rate. 

We compare the results of our simulations to the known 
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Mbh ~ & relation at z = (Eq. [5]) in Figure [7] and to 
the observed Mbh — Mtot relation at z = 0.2 (Eq. [5]) in 
Figure [8] We have arranged these plots such that the 
most "pessimistic" combinations of models — constant 
initial mass, distance and velocity tests for mergers, and 
new halo and merger checks evenly spaced in log a — are 
located in the lower-left portions of the plots, while the 
most "optimistic" scenarios are in the upper right. Our 
choice of model can greatly affect a number of aspects 
of SMBH relations, including the scaling of SMBH mass 
with a and M to t, the maximum mass of an SMBH, and 
the amount of scatter in the produced relations. 

Models with instant SMBH mergers produce relations 
with somewhat lower scatter, especially at higher masses, 
and generated higher maximum black hole mass than 
models with distance and velocity checks for mergers. 
For "prox" models, larger halos may contain several un- 
merged SMBHs, and spurious halo mergers do not lead 
to mergers of the black holes. These cause even the most 
massive halos to host relatively small black holes. Ad- 
ditionally, the difference in the total merger rate can be 
very dramatic, as shown in Figure^ The merger rate is 
defined as 



d 2 N 



AN 



dzdt AzAV 



47rc(l + zfd\{z), 



(8) 



where AN is the number of mergers in the redshift in- 
terval Az, AV is our simulated volume, c?a( z ) is the 
angular diameter distance, and c is the speed of light. 
The difference between models is especially significant 
at low redshift, where "prox" models can reduce the 
peak merger rate by a factor of two. The difference is 
negligible at high redshift, since the merger rate here is 
driven mostly by collisions of smaller halos, and the dif- 
ferences between "prox"- and "halo" -based merging are 
the smallest. For all cases, our measured m erger rate is 
less th a n rates found by m erger trees fe.g.. iMicic et al.l 
(|2008fk iMenou et all (|2001l )) since those works can in- 
clude SMBH masses below those which we can resolve. 

The frequency of merger checks can have a moderate 
impact on the slope of the final relations and a signifi- 
cant impact on the maximum black hole mass. However, 
changing this portion of the model does not significantly 
alter the scatter for models with constant initial mass. 
Checking for mergers every time step raises the maxi- 
mum SMBH mass by roughly half an order of magnitude 
and allows smaller halos (< 10 12 Mq) to host black holes 
fitting the observed relations. We can explain these re- 
sults by studying the merger rate in Figure |H1 Here, we 
see that merger checks at every time step can increase 
the merger rate by an order of magnitude relative to 
merger checks evenly spaced in log a. This is largely due 
to an increased number of seed black holes: by check- 
ing at every time step, we may capture halos as soon 
as they become resolvable, and before they encounter 
their first merger event. Indeed, we seed roughly twice 
as many SMBHs in the "dt" models relative to the "log" 
models, even though they both end up with roughly the 
same number of black holes at z = 0. Also, we capture 
more mergers at late times: in "log" models, we may skip 
the formation of small halos and their merging onto an 
already-formed larger neighbor, missing the SMBH mass 
associated with the smaller halo. This occurs regardless 
of seeding method. Spurious re-seeding also contributes 



somewhat, but this does not dominate because we see the 
same relative numbers of seeded and final black holes in 
both the "prox" and "halo" models. 

Using an initial seed based on early Mbh — cr relations 
produces broad scatter in the final relations. This is espe- 
cially evident in the "m-s, halo, log" combination, in which 
there appear to be two distinct populations of SMBHs: 
one population along the observed relation and another 
at lower mass. This behavior is due to the fact that 
halos are seeded in two scenarios: when the halo first be- 
comes resolvable in the simulation, and if the halo merges 
with another halo, loses its black hole, and later sepa- 
rates. Figure [TU] illustrates this by showing the initial 
seed mass as a function of redshift. Note that seed black 
holes are never larger than ~ 10~ 5 of the host halo mass. 
Small halos continue to appear throughout the evolution 
of the simulation and are seeded with ~ 10 4 Mq black 
holes. However, occasionally a larger halo loses its black 
hole and must be re-seeded with a correspondingly larger 
SMBH. Thus we may be left with two populations: halos 
with their original SMBHs that evolve to relations sim- 
ilar to halos with constant initial mass, and halos that 
are re-seeded at lower redshifts with black hole masses 
closer to the observed relation. This distinction is largely 
eliminated by applying the merger tests at every time 
step and merging black holes instantly on halo mergers. 
This results both in higher mass due to increased merg- 
ers and in increased low-redshift re-seeding of high-mass 
halos. This implies that re-seeding selects a few high- 
mass halos and places them on the observed relation at 
late times. Indeed, we find that all of the SMBHs with 
Mbh > 10 9 Mq and those black holes that lie at or above 
the observed relations (roughly half of the population 
at intermediate masses) are the product of low-redshift 
re-seeding. However, the majority of black holes evolve 
without re-seeding at late times, so that the differences 
between "m-s,dt" and "m-s, log" models at intermediate 
masses are largely due to the increased merger rate pro- 
moting the masses of all black holes. Note that the earlier 
discovery of halos by checking at every time step does not 
play a significant role here, since halos discovered later 
are simply initialized with larger black holes (which was 
one of the objects of this model). 

We may also understand these results by examin- 
ing the cosmic SM BH mass dens i ty, as in Figure QTJ 
Note that following Shanka r et al.l (j2004f ) we only count 
black holes that are matched to a halo and have masses 
10 6 < Mbh/ Mq < 5 x 10 9 when computing the mass 
density. All the SMBH densities produced by our mod- 
els are well below the observed z = density of 4.3 x 
10 5 M Mpc~ 3 (jShankar et al J 120041 ). Despite the fact 
that the smallest halos in the "m-s" runs started with 
10 4 Mq seed black holes, re-seeding and merging at low 
redshift causes the largest black holes to reach ~ 10 9 Mq 
and the cosmic mass density to reach nearly 1/10 of the 
observed value. 

Even though the re-seeding artificially creates all the 
SMBHs larger than 10 9 Mq, and half the black holes 
above 10 Mq, the majority of black holes are not re- 
seeded and have masses lower than that of the relation. 
Thus, even though we can reproduce the m-sigma re- 
lation at high masses (albeit with high scatter) there 
are not enough black holes with sufficient mass to reach 
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Figure 7. Comparison of models (shown as points) against the observed — a relation at z = (solid line). Dashed lines indicate la- 
uncertainty bounds in the observed relation. 



the observed cosmic mass density. Similarly, the higher 
merger rates produced by instant merging models pro- 
mote more black holes above the minimum mass thresh- 
old for inclusion in the density calculation. The fre- 
quency of merger checks has a dramatic impact on the 
z = density: here there is up to a factor of five dif- 
ference between models. This is largely due to the in- 
creased rate of discovering early halos in models with 
constant initial mass, thereby adding more SMBH mass 
to the simulation at early times. However, the extra mass 
added when checking "m-s" models at every time step is 
largely due to increased re-seeding. 
We quantify our comparisons by fitting our resulting 
— cr and Mbh — M to t relations to a straight line in 
log space. Since the "prox" models are obviously poor 
fits to power laws, we will not include them. To mimic 
observations, we only include black holes with masses 
greater than 10 6 Mq. Tables [3] and 0] describe the differ- 
ences between fits to our models and the observed rela- 
tions. We assume the halo mass uncertainty of Eq. (O; 
however, since it is difficult to quantify all the uncertain- 
ties for a generated in cosmological simulations, we will 
ignore them in computing fits. While this is an admit- 
tedly crude procedure, it does give us some estimate of 
the ability to distinguish these models from observations 
and from each other. We define the error in the tables 
to be the difference between the result of the model in 



Table 3 

Best-fit Slope and Normalization with lcr Uncertainty, Compared 
to Observed Afbh — CT Relation at z = (Eq. [6} . 



Model 


Slope 


Error 


Normalization 


Error 


con, halo, log 


2.2 ± 0.1 


-4.4 


6.9 ± 0.3 


-3.2 


con, halo, dt 


1.9 ± 0.1 


-5.1 


7.3 ± 0.4 


-2.0 


m-s, halo, log 


4.6 ± 0.3 


0.6 


8.0 ± 1.0 


-0.2 


m-s, halo, dt 


4.1 ± 0.2 


-0.3 


8.1 ± 0.5 


-0.1 



Note. — Error is defined in the text. 



Table 4 

Best-fit Slope and Normalization with lcr Uncertainty, Compared 
to Observed Afbh — Mtot Relation at z = 0.2 (Eq. 0. 



Model 


Slope 


Error 


Normalization 


Error 


con, halo, log 


0.89 ± 0.01 


-2.0 


6.8 ± 0.2 


-4.9 


con, halo, dt 


0.87 ± 0.01 


-2.1 


7.2 ± 0.2 


-3.4 


m-s, halo, log 


2.1 ± 0.2 


1.0 


7.0 ± 1.0 


-1.0 


m-s, halo, dt 


1.70 ± 0.05 


0.4 


7.9 ± 0.4 


-0.6 



Note. — Error is defined in the text. 



the simulation and the observed quantity divided by the 
sum of their respective la uncertainties. 

We see that, in general, models with constant initial 
mass are indistinguishable from each other, while these 
models are, as a group, significantly different from mod- 
els with varying initial mass. None of the constant-mass 




Figure 8. Comparison of models (shown as points) against the observed M^h — Mtot relation at z = 0.2 (solid line). Uncertainties are 
not shown but are included in later analysis. Dashed lines indicate la uncertainty bounds in the observed relation. 
0.12 



0.1 



0.06 



0.04 







halo.dt 






prox.dt ■ 






halo, log 






prox.log 









1 2 3 4 5 6 
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Figure 9. Merger rate of all SMBHs in the simulation volume as 
a function of redshift. We only show models with constant initial 
mass. Models with instant halo-based merging arc shown with 
solid lines, and models with merging tests are shown with dotted 
lines. Bold lines describe models with merger checks at every time 
step; thin lines show models with infrequent checks. 

models produce enough black hole mass to match the ob- 
served relations, as seen in the low normalization values, 
and they do not produce enough high-mass SMBHs in 
the largest halos, as seen in the low slope. 
In contrast, seeding halos with black holes with masses 




4 6 
Redshift 

Figure 10. Initial SMBH mass for seed black holes as a function 
of redshift for the "m-s,halo,dt" run. 

from an early Mbh — o relation matches observations well 
at low redshift. While these models fit the observed re- 
lations well, the total mass density of SMBHs is still well 
below the observed value. 

Surprisingly, we see marginally steeper slopes in the 
"log" models than in the "dt" models. Even though "dt" 
models have higher merger rates and produce more mas- 
sive black holes, they also generate a large population of 
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halo.dt 




Redshift 

Figure 11. Cosmic density of SMBHs with masses 10 6 M Q < 
M bh < 5 x 10 9 M© as a function of redshift. The top panel shows 
all models with constant initial seed mass, while the bottom panel 
shows models with seed masses based on the M^h — cr relation. 
Models with instant merging arc shown with bold lines; those with 
distance and velocity tests for merging are shown with thin lines. 
Finally, solid lines describe models with merger checks at every 
time step, and dashed lines indicate models with checks evenly 
spaced in log a. 

moderate-mass black holes hosted in small halos. This 
skews the fits to flatter slopes despite the higher max- 
imum black hole mass. However, these differences are 
not statistically significant, except for the constant ini- 
tial mass combined with instant merging scenarios, where 
the low scatter produces small uncertainties. 

Generally, the difference between a particular model 
and the Mbh — o relation is matched by the difference 
between that same model and the Mbh — M to t relation. 
However, the models are more indistinguishable from the 
Mbh — M tot relation, mostly due to the larger uncertain- 
ties in the observed relation. 

5. CONCLUSION 

We have developed a new fast, parallel halo finder for 
inclusion in cosmological simulations with the simulation 
code FLASH. Using SO halo finding techniques, we are 
able to produce halo catalogs in good agreement with 
traditional post-processing halo finders. Since our halo 
finder is designed to be fast, we are able to perform halo 
finding operations at every time step in the simulation, 
allowing us to perform a detailed analysis of SMBH sub- 
grid models. 

While merging alone cannot generate enough total 
mass in SMBHs to match the observed cosmic mass den- 
sity or generate high enough maximum black hole mass 
in the largest halos to match the observed Mbh — c and 
Mbh — Mtot relations, it can play a large role in devel- 
oping the slope of the relations, especially at intermedi- 
ate mass ranges. Thus, merging should not be totally 
discounted in considering the processes that provide the 
correlations between black hole mass and bulge, galaxy, 
and halo properties. However, since none of our consid- 
ered models can account for the observed cosmic mass 
density, there is still a significant role for accretion and 
feedback processes in the evolution of SMBHs. Also, 
since the choice of models can greatly affect the cosmic 
SMBH mass density, accretion and feedback models must 
be chosen carefully to match observations. 



The choice of subgrid models can dramatically impact 
the merging rate of black holes. Since the merger rate 
has a large influence on the performance of upcoming 
gravitational wave detectors, halo finding operations in 
simulations should be done as frequently as possible in 
order to accurately capture this rate. We do not believe 
that the inclusion of gas will significantly affect these 
relative differences, since they are largely driven by the 
ability to find more black holes at early times. However, 
the choice of merging test does not greatly affect the 
predicted rate, except at low redshift. 

While we have bracketed the possible subgrid mod- 
els with "optimistic" and "pessimistic" scenarios, models 
best matching insights from theory and observations are 
usually in between those extremes. While seeding black 
holes with a uniform initial mass for black holes may well 
model high-rcdshift behavior, it is not clear that this is 
a useful strategy for re-seeding low-redshift halos. Most 
re-seeding is certainly an artifact of the halo finder and 
the lack of an in-code merger tree to determine when ha- 
los have truly merged. However, there are some plausible 
scenarios where re-seeding may be needed: for example, 
when three-body or gas interactions strip an SMBH from 
a central galaxy. Also, when this approach is coupled 
with infrequent halo finding operations, it may deposit 
too little mass in the seed SMBHs. While instant merg- 
ing is too optimistic, our current lack of understanding 
of SMBH mergers implies that we cannot entirely specify 
this portion of the subgrid model, and we must rely on 
a bracketing procedure. 

Future examination of these subgrid models must be 
done in a cosmological simulation involving gas evolu- 
tion, accretion onto the black holes, and feedback from 
active galactic nuclei. However, since models of these 
processes carry with them their own assumptions and 
adjustable parameters, care must be taken to fully sep- 
arate the effects of black hole seeding and merging. It 
may be possible that due to the self-regulating nature of 
feedback that the differences among these models may 
disappear; however, the differences in merger rate and 
peak black hole mass suggest significant variances may 
remain. Only once all aspects of these subgrid models 
are analyzed, understood, and compared to our observa- 
tional understanding can we confidently combine them 
into an integrated model. 

ACKNOWLEDGMENTS 

The authors acknowledge support under a Presidential 
Early Career Award from the U.S. Department of En- 
ergy, Lawrence Livermore National Laboratory (contract 
B532720). Additional support was provided by a DOE 
Computational Science Graduate Fellowship (DE-FG02- 
97ER25308) and the National Center for Supercomput- 
ing Applications. The software used in this work was 
in part developed by the DOE-supported ASC/ Alliance 
Center for Astrophysical Thermonuclear Flashes at the 
University of Chicago. This research used resources of 
the National Center for Computational Sciences at Oak 
Ridge National Laboratory, which is supported by the 
Office of Science of the US Department of Energy under 
contract no. DE-AC05-00OR22725. The authors thank 
Brian O'Shea, Tiziana Di Mattco, and many others for 
useful conversations at the tenth Great Lakes Cosmology 
Workshop. 



11 




10 100 1000 10000 100 1000 

Number of Processors Number of Processors 

(a) Strong scaling of pSO halo finder at various problem sizes as a (b) Weak scaling at smallest number of processors for each problem 
fraction of total wall clock time in a single FLASH time step. size as a fraction of total wall clock time in a single FLASH time 

step. The dotted line indicates weak scaling after adjustments have 
been made to account for non-linear scaling of halo counts with 
problem size. 

Figure 12. Strong and weak scaling of the pSO halo finder. 



APPENDIX 

PERFORMANCE AND PARALLEL SCALABILITY 

In order to minimize communication among processors, we divide potential halos into two lists. For every potential 
halo, we compute the overdensity at the largest possible on-processor radius. If the next search radius is larger than 
this, the halo is added to a list that must be communicated. Thus, all on-processor halos are processed concurrently 
with no communication. Since the volume of communication is very small, we send all other halos to every processor. 
However, since this approach is not scalable to very large halo catalogs or numbers of processors, we switch to a 
buffered communication pattern if we cannot allocate a global halo catalog. In the buffered approach, halo lists are 
only communicated to the nearest processors one at a time. The list is communicated until all blocks within the 
volume of the processor's halos have been searched. While slower than the all-to-all approach, the amount of storage 
per processor required stays fixed as the number of processors and number of halos grow. In either case, communication 
continue s until all halos are fully searched. For this study, we will use the all-to-all approach. 



Figure 12(a) shows the strong scaling of the pSO halo finder as a fraction of the total wall clock time in a single 
FLASH time step. We show three different uniform problem sizes: 256 3 , 512 3 , and 1024 3 particles and zones. We 
restarted each simulation at a representative redshift, z = 0.25, and ran for five time steps. These times do not include 
optional portions of the halo finder routine, such as writing the halo catalog to disk or tagging particles within halos. 
We performed these calculations on jaguar, a Cray XT5 system at Oak Ridge National Laboratory. Jaguar consists 
of 16,688 dual six-core AMD Opteron nodes with 16 GB of memory per node and has a peak performance of 2.332 
petaflops. Our approach offers good strong scaling behavior: we are able to beat or match the scaling performance 
of FLASH at all problem sizes. At larger core counts, FLASH has difficulty scaling the Poisson solver, whereas the 
halo finder maintains good scalability. However, we can infer that we have poor weak scaling: the halo finding steps 
require ever larger wall clock time as the problem size grows. 

The poor weak scaling is due to several factors. First, since these are uniform grid calculations, at low redshift the 
particle distribution among processors becomes highly unbalanced. FLASH is block based and uses a Morton curve 
the distribute these blocks among the processors. Thus, while each processor has roughly the same number of blocks, 
those blocks in highly dense regions will contain many more particles than those in voids. Since our halo finder scans 
through particles, there is a lack of concurrency due to this imbalance. At small problem sizes, this is not an issue, but 
with 1024 3 particles all processors must wait for the few heavily loaded processors to complete searching. There is an 
imbalance in halos as well, since highly overdense regions contain more halos than underdense regions. Note that this 
is a bigger problem for box sizes smaller than 200 Mpc. We found that the 256 3 run placed a maximum of 15 halos on 
a single processor, the 512 3 r un had up to 33 halos on a single processor, and the 1024 3 run mapped up to 80 halos 



on a single processor. Figure 12(b) shows the effects of the imbalanced load due to an increased maximum number 
of halos on the weak scaling behavior. We compute the effects by maintaining no more than 15 halos per processor 
no matter the problem size. We performed this test by ignoring halos past this limit. We see that halo imbalance 
explains a significant part of the weak scaling results. 

Second, as we increase the resolution, we greatly increase the number of resolvable halos. Since smaller halos are 
much more numerous than larger ones (see Figure [SJ , for every doubling of the resolution we get roughly ten times 
the number of halos. However, the number of particles and zones only increases by a factor of eight. This means that 
even if we devote eight times the number of processors to the problem, we will spend a larger fraction of a time step 
finding all the halos. 
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Finally, as we increase the resolution, we require more steps in the binary search procedure to meet the stopping 
criterion. At 256 3 zones, we require roughly 35 steps to identify the halos. At 512 3 zones, we require 45 steps on 
average. Finally, for 1024 3 zones, we average 55 steps. This increases the serial runtime of the halo finder, regardless 
of the degree of concurrency. 

These calculations were performed with only dark matter. As we perform runs with more included physics, such as 
hydrodynamics and radiative cooling, the fraction of time spent in halo finding will decrease, even for high-resolution 
runs. Also, we expect AMR runs to alleviate the problems of poor particle load balancing, since in cosmological 
simulations we typically refine on overdense regions. We are currently preparing a manuscript detailing the performance 
of the halo finder with different box sizes, combinations of resolutions of particles and zones, additional physics, and 
different AMR refinement schemes. 
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